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I. INTRODUCTION 


With the formation of the NASA Office of Exploration and the 
July 20, 1989 speech by President George Bush initiating a space 
exploration planning study, more attention has been focused on the 
planning and technologies required to establish a base on the lunar 
surface. Many different technologies have been considered to 
support such an endeavor with a great deal of attention being given 
to transportation systems. Many past studies of various lunar base 
strategies have led to the conclusion that a major cost driver of a 
base on the moon would be the transportation of the infrastructure 
from Earth to the lunar surface. One element of the transportation 
system is the transfer of crew and cargo between low-Earth orbit 
and lunar orbit and the subsequent return of such payloads. Many 
technologies have been applied to the space transportation in order 
to reduce the program costs. One such technology that may be 
applied to orbital transfer is that of constant, low-thrust propulsion 
systems. These systems, generally more efficient than conventional 
chemical propulsion systems tend to reduce the fuel requirement by 
several orders of magnitude and thus reduce the Earth to orbit 
launch requirements which results in a cost savings to a lunar base 
program. However, it must be noted that low-thrust propulsion 
systems are not applicable to crew transfer as the longer trip times 
would expose the crew to the hazards of the Van Allen radiation 
belts and long-term degradation of human physiological systems. 
Therefore, most studies addressing the low-thrust transfer of 
payloads from low-Earth orbit to lunar orbit have been primarily in 
unmanned cargo vehicles. This is the primary focus of this study. 

The analysis contained herein focuses on both the definition of 
the vehicle and its supporting systems in conjunction with the 
guidance and control algorithms required to fly the spacecraft from 
its departure orbit to its destination orbit. The vehicle definition and 
systems characterization portion of this study concentrated on 
collecting state-of-the-art information from such institutions as 
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NASA Lewis Research Center, Los Alamos National Laboratory, the 
Auburn Space Power Institute, United States Air Force Space 
Division, and the Strategic Defense Initiative Organization in order to 
develop sizing algorithms to characterize the most state-of-the-art 
spacecraft for assessment of its performance. The trajectory analysis 
part of this study is the primary state-of-the-art advancement of 
this study as it is a three dimensional analysis formulated in the 
restricted three-body problem with perturbations. 

In the past there has been a number of studies conducted on 
the low-thrust spacecraft, but most of these studies have either been 
dedicated to system characterization or trajectory analysis with very 
few being a combination of the two. Most of the trajectory analyses 
performed have been done in in-plane without out-of-plane 
thrusting required to change inclination of line of ascending node. 
This leads to a good initial characterization of a vehicle, but hardly 
lends itself to what is required for a complete, robust conceptual 
design and mission plan. Other past studies have also concentrated 
on a centric trajectory analysis. That is, a problem formulation in a 
single gravitational field with a spacecraft in orbit and spiraling 
outward (or inward) until escape (or capture) and then changing to a 
different centric coordinate systems for continued analysis. It should 
be noted that the Earth-Moon systems does not lend itself to this 
type of analysis as the Moon is in close proximity to the Earth's 
sphere of influence. Therefore, the problem formulation was 
developed in the restricted three-body problem to facilitate a more 
accurate simulation of the spacecraft's operation. 

Since the trajectory analysis technique chosen here is based on 
a numerical integration, the vehicle must be initially sized before the 
trajectory analysis can begin. Therefore, an initial estimation model 
for the spacecraft size was developed which then transfers data to 
the trajectory analysis model. After the trajectory analysis is 
completed, the spacecraft is re-sized to fit the computed trajectory. 
This iteration scheme can be repeated until a suitable level of fidelity 
is reached. 
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II. PROBLEM FORMULATION 


To facilitate the development of the low-thrust cislunar 
spacecraft trajectory determination and sizing model it was 
important to segment the total problem into discrete problems that 
could be more easily addressed. From a top-level perspective, the 
problem was divided into two functional areas: system 
characterization, and trajectory generation and guidance. These two 
functional areas were further divided to better enable a solution to 
determined. 

In the case of the system characterization portion of the model 
it was necessary to determine how multiple technologies for various 
components were to be integrated into a complete system. Also of 
concern was the determination of the amount of propellant required 
by the vehicle. The propellant requirement is needed to enable the 
determination of an appropriate initial mass of the spacecraft in the 
departure orbit. The value of the spacecraft's initial mass is used to 
determine the initial acceleration of the spacecraft and the 
subsequent acceleration throughout the trajectory generation. 

For the trajectory generation and guidance segment of the 
model, two reference systems were chosen to allow the greatest ease 
of numerical integration, conceptual understanding, and control for 
various portions of the spacecraft’s flight. The first reference system 
is a non-dimensionalized geocentric frame in which six orbital 
elements describe the shape and nature of the spacecraft's orbit. 

This is the equinoctial coordinate system with modifications. The 
second system is a rotating, non-dimensionalized, right handed 
coordinate system with the origin at the Earth-Moon system center 
of mass. This coordinate system is known as the restricted three- 
body system. Both reference systems are explained in later sections. 

2.1 Vehicle System Characterization 

The characterization of a low-thrust orbital transfer vehicle 
(OTV) for Earth-Moon transfers was performed by dividing the 
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spacecraft into its major systems. These are power, propulsion, 
thermal control, propellant tankage, and structural mass. Each of 
these systems is defined in terms of interactions between the other 
systems. For example, the propulsion systems needs power, thermal 
control, propellant, and structure, the power system needs thermal 
control, and structure, the thermal control system needs power, and 
structure, etc.... 

The payload mass is the primary driver when characterizing 
the OTV. The user can input any desired size for the payload and the 
system sizing model will characterize a spacecraft based on the 
payload mass. The propulsion system is chosen by the user to define 
the performance of the system. The propulsion system performance 
characteristics, such as specific impulse, and power per engine, along 
with the technology choice of the propulsion system are required 
user inputs. These inputs allow the model to adequately define the 
propulsion system of the OTV. The propulsion system then has 
certain requirements that the other systems must provide in order to 
have a functioning spacecraft. The technology choice of the power 
system is another user input that completes the necessary user 
characterization of the OTV. The power system obtains the 
information regarding its electrical requirements primarily from the 
propulsion system. The thermal control system obtains the required 
heat rejection data from the propulsion and power systems and the 
structural requirements are developed from the mass and volume 
requirements of the propulsion, power, and thermal systems. 

Sizing the subsystems required for an OTV is a formidable task 
due to the large amount of data needed to cover all of the possible 
system choices. In order to implement a computer-based vehicle 
systems sizing model, it was determined that the use of parametric 
equations would be more effective than using a large design 
database to generate system sizing parameters. The use of 
parametric equations results in a robust design model. For example, 
the parametric equations used in the power system sizing subroutine 
were derived by determining mathematical relationships between 
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such system parameters as output power and system mass. Using 
these relationships, the power system mass can be calculated for any 
desired output power within the range of the power system 
technology. All of the system sizing subroutines use parametric 
equations, derived from data in current technical publications, to 
generate the system mass estimates. 

2.2 Trajectory Generation 

In previous studies on the development of trajectories for low- 
thrust cislunar OTVs, little attention has been directed toward the 
robust algorithms required for three-dimensional guidance and 
control of the spacecraft. The guidance and control of the spacecraft 
and the determination of the appropriate trajectory are highly 
coupled. The guidance, control, and trajectory determination are 
closely related problems which by necessity must be treated with 
equal importance 1 . A major problem in the design of low-thrust 
OTVs and their associated trajectories is the lack of an end-to-end 
computer simulation for the spacecraft trajectory between Earth and 
lunar orbit. The physical capability of a low-thrust spacecraft to 
travel between the Earth and the Moon is known 2 . The current 
question is how the vehicle will behave at the proposed thrust level 
and how it will be guided on its trajectory. To adequately model the 
dynamic forces on the spacecraft in its travel in the Earth-Moon 
system, two problem formulations were used to generate the 
trajectory. The equinoctial formulation of the equations of motion 
was used for Earth escape, and the restricted three-body formulation 
was used for the for the midcourse and capture phases. 

The standard formulation of a spacecraft’s orbit is through the 
use of classical elements. These correspond to the semi-major axis of 
the orbit (a), eccentricity (e), inclination (i), mean anomaly (M), 
longitude of periapsis (to), and longitude of the ascending node (f2). 
These classical elements are usually defined in terms of a geocentric 
equatorial coordinate system with the positive X-axis parallel to a 
vector from the sun to the Earth at vernal equinox, and the Z-axis 
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through the north pole of the Earth (see Figure 1). These six 
elements are traditionally used to describe the orbit of a spacecraft. 

In the classical formulation of the two-body problem, the 
angles of mean anomaly (M), eccentric anomaly (E), and true 
anomaly (v), are well known and much used. Low-thrust spacecraft 
trajectories will always have portions of their trajectories where the 
orbit is circular or the eccentricity is very small («1). Also, it is 
possible to have small to zero inclination orbits when transferring to 
geosynchronous orbits or passing through the equatorial plane. The 
three angles, M, E, v, are not well suited for small eccentricity or 
circular orbits because they are measured from the perigee or 
periapsis point of the orbit. When a orbit is circular the periapsis 
point of the orbit is not defined. Additionally, © and Q are ill- 

defined for zero inclination orbits. This causes singularities during 
numerical integration of orbits with classical orbital elements. For 
these reasons, integration of low-thrust trajectories in classical 
orbital elements is not feasible or desirable. 

2.2.1 Equinoctial Formulation 

An alternate set of orbital elements can eliminate the 
singularities experienced by classical orbital elements. These 
elements are known an equinoctial elements and are defined in 
terms of classical elements, 

a = a, 

h = e sin (© + Q), 
k = e cos (co + Q), 

X = M + © + £2, 
p = tan i/2 sin Q, and 
q = tan i/2 cos Q. 

The equinoctial element formulation, p a = (a, h, k, X, p, q), has 
no singularities at e=0 or i= 0. The semi-major axis, a, remains the 
same while h, k, p, and q are just convenient mathematical 
formulations, and X is the sum of M, E, and v and is referred to as the 
mean longitude. 
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The equinoctial coordinate frame is defined by unit vectors f, g, 
and w illustrated in Figure 2 and defined as follows. 

1 - p 2 + q 2 

f = L 


1 + p 2 + q 2 


2pq 

-2p 


g = 


1 + p 2 + q 2 


2pq 

1 + p 2 - q 2 
2q 


w 


1 

1 + p 2 + q 2 


2p 

-2q 



Where each column vector in brackets contains the X, Y, and Z 

A 

component of the unit vectors f, g, and w. 

The trajectory of a low-thrust OTV can be integrated in this 
formulation through the use of a method known as variation of 


parameters. A complete derivation of this method is included in 
Appendix A. The equinoctial elements are integrated using the 
variation of parameters equation written in equinoctial form: 

P« = %-F, 

a; {1) 

where the right hand side of { 1 } is the dot product of the partials 
with respect to velocity jvith the rectangular components of the 
perturbing acceleration F. For the low-thrust spacecraft, the 
perturbing acceleration vector, F, consists of the acceleration due to 
the spacecraft s propulsion system, the perturbing acceleration of the 
Moon, the gravitational pull due to the oblateness of the Earth, and 
the acceleration on the spacecraft due to atmospheric drag. 

The equinoctial elements change slowly during integration and 
are smoothly varying functions. This allows the numerical integrator 
to take large step sizes, thereby speeding calculation of the 
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Figure 2 - Equinoctial Coordinate Frame 
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trajectory. For both the increase in integration speed and the 
elimination of singularities, the equinoctial formulation is used to 
generate the trajectory during Earth escape. 

2.2.2 Restricted Three-Body Formulation 

To adequately understand the dynamics of motion of the low- 
thrust spacecraft, the gravitational effects of the Earth and the Moon 
on the spacecraft must be included for the full duration of the 
trajectory. The thrusting acceleration of low-thrust OTVs in high 
Earth orbit is the same order of magnitude as the perturbing force 
due to the Moon. To model the Earth-Moon system with the 
necessary accuracy and achieve computational efficiency, the 
restricted three-body formulation of the dynamical equations is 
utilized as the governing equations of motion. 

The problem of three bodies was first formulated in 1772 by 
Lagrange. Further studies by Poincare, Laplace, Hill and Szebehely 
have resulted in a detailed treatment of the problem and a general 
understanding of the interactions between the two primary 
gravitational fields. Various formulations are available to represent 
the problem of three bodies. Discussions with Victor Szebehely at 
The University of Texas at Austin led the authors to use the non- 
dimensional formulation of the restricted problem of three-bodies 
for the cislunar transfer and capture portion of the trajectory 
generation 3 . This formulation has several advantages over the 
general three-body problem. The order of equations to be integrated 
are reduced from 18th order for the general problem of three bodies 
to 6th order for the restricted problem of three-bodies. This 
reduction in order dramatically decreases integration time. 
Additionally, the equations of motion are non-dimensionalized by the 
Earth-Moon distance and the angular rate of the Earth-Moon system 
about the systems center of mass (barycenter). The Earth-Moon 
system in restricted three-body formulation is shown in Figure 3. 

Many realistic orbital cases of interest permit treatment as 
restricted three-body situations. A case of particular interest is that 
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Figure 3 - Restricted *] 
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of a spacecraft moving in the Earth-Moon system. Certain 
assumptions can be made about the nature of the Earth-Moon system 
that permit a simplified formulation at a slight loss of accuracy. It is 
assumed that the motion of both the Earth and Moon is circular and 
coplanar about their barycenter. The Earth-Moon line is the x-axis of 
a rotating coordinate system and the z-axis is parallel to the angular 
momentum vector of the Earth-Moon plane. The y-axis is 
perpendicular to both the x-axis and the z-axis. The position of the 
Earth and the Moon are constant in the restricted three-body 
formulation. The spacecraft, at point P, is assumed to have negligible 
mass and to have no impact on the motion of the Earth or the Moon. 
The motion of the spacecraft is governed by the relative gravitational 
attraction of the Earth and the Moon rotating about the barycenter. 
The equations of motion for the spacecraft in the restricted three- 
body non-dimensionalized rotating coordinate system are: 



x - 2 y = £2 X , 

{2} 


y + 2x = C2y , and 

{3} 


i= n z , 

{4} 


2 r \ r 2 , 

l/ 

{5} 

r \ = 

(x - /i) 2 + y 2 + z 2 ] 2 , and 

/ 


r 2 = 

[(x + 1 - /i) 2 + y 2 + z 2 ] 1 . 



A derivation of the equations of motion for the restricted three-body 
problem can be found in Victor Szebehely's book. Theory of Orbits A . 

In the restricted three-body formulation of the equations of 
motion the Keplerian (potential and kinetic) energy of the spacecraft 
is not conserved. This is due to the assumptions regarding the 
motion of the Earth and Moon and the spacecraft's effect upon their 
motion. However, the sum of the angular momentum, velocity in the 
rotating coordinate system, and potential energy of the spacecraft is 
conserved. This result is determined through the derivation of the 
equation known as the Jacobian integral. The integral and the 
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constant, C, it produces are named after mathematician Karl Gustav 
Jacobi who first formulated this integral in 1836. This integral can 
be derived from equations {2}, {3}, {4}, and {5}. Initially, equations 
{2}, {3}, and {4} are multiplied byx ,y , and z respectively. This 
yields, 

x x - 2y x = Q x x , 

y y + 2x y = Gy y , and 
z i- z • 

These three equations are added together to form, 

x x + yy+ z z = Cl x x + Q y y + £2 Z i , 

which can be integrated and rearranged to find the integral, 

C = 2n-(x 2 + y 2 +z 2 ). {6} 

This constant, C, can be determined for any set of position and 
velocity conditions in the restricted three-body problem. If there is 
no force acting on the spacecraft other than the Earth and Moon the 
value of C, the Jacobian constant, will remain constant. So for a 
nonperturbed orbit about the Earth the value of C will remain the 
same. For various arbitrary combinations of position and velocity 
differing values of C will be found. Conversely, if C is initially 
determined to be a particular value from equation {6}, say the value 
that corresponds to the position and velocity of a spacecraft in orbit 
about the Earth, there will be many combinations of position and 
velocity the will give the same value of C. The potential, Q, shown in 
equation {5} is dependent solely upon the position of the spacecraft 
in the restricted three-body system. If x , y , and z were each set 
equal to zero in equation {6} and a value of C was previously 
determined from a spacecraft's position and velocity, then equation 
{5} will give the x , y, and z coordinates where the spacecraft would 
have zero velocity. These coordinates form a surface in the three- 
dimensional coordinate system of the restricted three-body 
formulation. This surface bounds the locus of points where the 
spacecraft can have motion given the value of its Jacobian constant, 
C, in the Earth-Moon system. On this surface a spacecraft with a 
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given C will have zero velocity. Only on the inside of the curve will 
the spacecraft have any velocity, thus restricting the motion of the 
vehicle to the inside. Figure 4 shows a series of cross sections of zero 
velocity surfaces as curves in the Earth-Moon system's x-y plane 5 . 
The value of the spacecraft's Jacobian constant is used during the 
midcourse targeting phase of the trajectory generation. It is used as 
an indicator of the spacecraft's ability to achieve the desired cislunar 
transfer. 

* Battin, R. H., Miller, J. S., “Trajectories and Guidance Theory for a Continuous 
Low-Thrust Lunar Reconnaissance Vehicle,” 6th Symposium on Ballistic 
Missile and Aerospace Technology, 1961. 

2 Hill, P.G., and Peterson, C.R., Mechanics and Thermodynamics of Propulsion . 
Addison-Wesley Publishing Company, Inc., Reading, Massachusetts, 1965, 
Chapter 10. 

2 Szebehely, Victor, Personal Communication, October 1988. 

4 Szebehely, Victor, Theory of Orbits . Academic Press, Inc., New York, 1967 

5 Kaplan, Marshall H„ Modem Spacecraft Dynamics and Control. John Wiley 
and Sons, New York, 1976. 
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Figure 4 - Zero Velocity Curves in the Earth-Moon System 
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Ill, PROGRAM METHODOLOGY 


3.1 Vehicle Systems Sizing 

The vehicle system sizing model contains parametric models of 
the major vehicle subsystems: power generation and distribution, 
propulsion, support and propellant tankage structure, and thermal 
control. In some cases multiple technologies are parameterized for 
each system with the technology choice left to the model user. The 
interrelationships between these systems are mapped out as a set of 
iterating functions that are scaled according to the users system 
technology and payload size choice. Figure 5 shows some of the 
interrelationships between the various systems. It is assumed that 
the OTV does not refuel at the Moon and that no propellant tanks are 
discarded. The total dry mass of the OTV is calculated by summing 
the system component masses and desired payload mass. 

A propellant estimation model has been developed to 
determine the appropriate amount of propellant necessary for a 
cislunar transfer. A functional relation between the spacecraft's 
thrust to weight ratio and the propellant fraction of the vehicle has 
been developed using a two-dimensional constant thrust lunar 
trajectory program developed by the Large Scale Programs Institute 
under a grant from the NASA Johnson Space Center. Using this 
functional relationship the propellant mass can be determined which 
leads to an initial estimate of the total system mass. The initial total 
system mass of the spacecraft, along with the propulsion parameters 
(mass flow rate and specific impulse (Isp)), is used as model drivers 
by the trajectory generation model. 

A description and overall program flow of the vehicle system 
sizing model is detailed in the next section. Following the program 
flow description, the specific formulation of the propellant estimation 
routine is discussed and assumptions concerning the development of 
the methodology are detailed. 
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SYMBOLS USED: 
m-1 = payload mass 
P-1 = payload power req. 
m() = mass of system 
P() = power requirements 
or power output of sys. 
T() = thermal requirements 
of system 


NOTE: System masses 
will be used to estimate 
mass of structure. 


INPUT: mdot, Isp, 
m-1, dV-imp, P-1 

mdot 

Isp 

PROPl 

SYSTEM 

1 

JLSION 
1 MODEL 

1 


P(payload) 


[m (prop, sys.)] 

P(prop.sys.) Estimated therm.cont.pow. 
requirement 


P(etc.) 


POWER 

SYSTEM MODEL 


T(pow.sys. 


rejection) 


[m(pow.sys.)j 


T(pow.sys. rejection) 


T(etc.) 


THERMAL CONTROLj 
SYSTEM 
MODEL 


-P(therm.cont-) 


[m(therm.cont.)j 


[m(etc.)j 


Sum of system masses 


Figure 5 - Basic OTV System Interrelationships 


17 




3.1.1 System Model Program Flow 

The vehicle system sizing model is segmented into input and 
systems sizing subroutines. The first subroutine run is the input 
subroutine, INPUTO, which checks whether the user wishes to run 
the vehicle system sizing model or specify his own overall vehicle 
characteristics. If the user wishes to specify the specific vehicle 
characteristics, the program prompts the user to enter the initial 
spacecraft mass, the total mass flow rate of the propulsion system, 
the specific impulse of the propulsion system, and the final mass of 
the vehicle. These parameters are transferred to the trajectory 
generation subroutines. 

If the user chooses to run the vehicle system sizing model, then 
the user has the choice of generating the spacecraft’s characteristics 
for a one-way or two-way trip. A one-way trip would consist of a 
Earth to Moon, or Moon to Earth, transfer. A two-way trip is an 
Earth-Moon-Earth transfer. At present, even though the system 
sizing model can generate characteristics for a vehicle supporting 
two-way transfers, the trajectory generation routine only supports 
one-way transfers. The characteristics of the spacecraft are output 
to a file, LOWTHRST.DAT for later evaluation. After these choices 
have been made, program control is passed to the subroutine, 
SYSMOD (SYStem MODel), which handles the calls to the various 
system sizing subroutines. The top level flow of data is shown in 
figure 6. 

SYSMOD begins by calling the input subroutine, INPUT1, which 
lets the user enter the mass of the payload to be delivered (and 
returned in the case of the two-way transfer). Next, SYSMOD calls 
the propulsion system sizing subroutine, PROPMOD (PROPulsion 
MODel), which allows the user to choose one of four types of electric 
propulsion systems. The four types of systems available include ion, 
magnetoplasmadynamic (MPD), arcjet, and a user specified system. 
There are three options for Ion propulsion: Xenon, Krypton, or Argon 
propellant. PROPMOD calls one of the four propulsion system models: 
ION, MPD, ARCJET, or USER1. The subroutines, ION, MPD, and ARCJET, 



Input Payload 
Mass 


Choose Propulsion System 


Arcjet 


MPD 


1 INPUT 
Isp, number of thrusters, 
and power input to each 
thruster 


Size propulsion system based on 
user inputs and parametric 
equations 


Choose Power System Technology 


Size power system and 
remaining vehicle systems 
using parametric equations 


OUTPUT 

mass flow rate, Isp, and 
vehicle "wet" mass to 
trajectory generation 
subroutines 


OUTPUT 

vehicle system mass breakdown and 
system characteristics for user 
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allow the user to specify the specific impulse, number of thrusters, 
and power input to each thruster. These input quantities are used in 
parametric equations to calculate the efficiency of the thruster 
system, the dry mass of the propulsion system (i.e. no propellant), 
and the total power requirements of the propulsion system. The 
subroutine, USER1, allows the user to specify the specific impulse, 
mass flow rate per thruster, power input to each thruster, number of 
thrusters, efficiency, and the dry mass of the propulsion system. 

Once the propulsion system has been specified, control of the 
program is returned to SYSMOD which then calls the power system 
sizing subroutine, PWRMOD. 

PWRMOD allows the user to choose one of six power generation 
and conversion systems 1 - 2 . The possible choices include: 

( 1 ) Liquid Metal Reactor using Rankine Cycle 
Conversion (1.5 kWe - 50 MWe), 

(2) Liquid Metal Reactor- NERVA Derivative using 
Closed Brayton Cycle Conversion 

(1.5 kWe - 50 MWe), 

(3) Solid Core Reactor using In-Core Thermionic 
Conversion (10 kWe - 50 MWe), 

(4) Liquid Metal Reactor using AMTEC Thermoelectric 
Conversion (1 kWe - 50 MWe), 

(5) NERVA Derivative Reactor using 
Magnetohydrodynamic Conversion 
(100 kWe - 100 MWe), and 

(6) SP-100 reactor with Thermionic conversion 
(100 kWe - 500 kWe). 

Once the power system is chosen, PWRMOD uses parametric 
equations to size the reactor, convertor, and control systems based on 
the power required. In addition, PWRMOD computes the efficiency of 
the power system and returns program control to SYSMOD to call the 
thermal control system sizing model, THERM. THERM computes the 
mass of the thermal control system based on parametric equations 
concerning the efficiency and power output of the power system. 
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The final systems sizing subroutine called by SYSMOD is the 
reaction control system sizing subroutine, RCS. RCS uses parametric 
equations to estimate the mass of a reaction control system and the 
propellant required to give the vehicle an initial estimate of a 100 
meter per second velocity change. The vehicle mass used for this 
estimate consists of the system masses already calculated, an 
estimate of the vehicle's structural mass, and a rough estimate of the 
propellant required for the trip. Once this system has been sized, the 
only remaining quantity that must be determined is the amount of 
propellant required for the trip. The subroutine that contains the 
propellant estimation algorithm is called PRPEST. 

3.1.2 Propellant Estimation 

The subroutine PRPEST (PRoPellant ESTimator) computes an 
estimate of the mass of the propellant required for a one-way trip. 
The inputs to the subroutine consist of the specific impulse, mass 
flow rate per thruster, the number of thrusters in the propulsion 
system, the vehicle initial mass (structural mass plus propulsion, 
power, thermal control, and RCS system masses), and the payload 
mass for the one-way flight. The subroutine uses the mass flow rate 
per thruster and number of thrusters to determine the total mass 
flow rate for the propulsion system. 

PRPEST runs an iterative loop to calculate an estimate of the 
propellant mass required for the cislunar flight. The vehicle's 
velocity is initialized to zero, the mass of the spacecraft is the 
vehicle's initial mass (including the payload mass), and the 
propellant mass is zero. The total velocity change, based on Apollo 
17 data, required to travel between the Earth and Moon is 9000 
meters per second 3 . PRPEST calculates an incremental Av to be 

given to the vehicle at each time step in the loop. The equation used 
to calculate this Av is: 

Av = mdot*g*Isp*dt/m {7} 

where, 

mdot = mass flow rate (kg/s). 
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g = constant of gravity (km.sec A 2), 

Isp = specific impulse of the engines (seconds), 
dt = time step (100 seconds), and 
m = instantaneous spacecraft mass (kg). 

The calculated incremental Av from equation 7 at each time step is 
added to the vehicle’s current velocity to obtain an updated velocity. 
The propellant mass is incremented by the total mass flow rate over 
the 100 second time increment. The propellant tank mass is 
incremented at each time step using a parametric equation that 
relates tank mass to propellant mass. The new propellant and tank 
masses are added to the vehicle's total mass to obtain a new value 
for the vehicle's instantaneous mass. These equations are repeated 
until the vehicle's velocity meets or exceeds 9000 m/s. 

In order to use PRPEST to estimate the mass of the propellant 
required for a two-way trip, SYSMOD first generates the total 
spacecraft mass for the “second leg” of the two-way transfer. To 
determine the propellant mass for the “first leg” of the transfer, the 
propellant mass and tankage is assumed to be part of the first leg’s 
payload mass. The propellant mass for the first leg is summed with 
the propellant mass for the second leg to obtain the total two-way 
propellant required. 

3.2 Trajectory Generation 

The trajectory determination and guidance of the cislunar low- 
thrust OTV is divided into three distinct phases (see Figure 7): 

Orbital Plane Alignment, Midcourse Targeting, and Capture and 
Circularization. Orbital Plane Alignment is concerned with the 
methodology and required guidance and control to drive the 
spacecraft from its initial orbital plane about the Earth into the plane 
of motion of the Moon about the Earth. The Midcourse Targeting 
phase of the trajectory generation deals with achieving the 
spacecraft position and velocity at Earth or Moon escape to achieve 
Earth-Moon transfer. The Capture and Circularization phase consists 
of the required controls to capture the spacecraft about the target 
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Figure 7 - Trajectory Generation Phases 




planet and to circularize the capture orbit at the desired final 
altitude above the surface of the target body. Each of these phases 
has a different guidance scheme to achieve the overall goal of 
generating trajectories between the Earth and Moon. 

3.2.1 Orbital Plane Alignment 

The first phase in any cislunar journey for an OTV is the escape 
from the initial parking orbit, whether about the Earth or the Moon. 
For low-thrust spacecraft to achieve escape, a long period of 
continuous thrusting is necessary. This results in a slowly increasing 
spiral trajectory from the initial orbit. It can be shown that a near- 
optimal thrust for a planar orbital transfer should be directed along 
the velocity vector of the spacecraft for the majority of the 
trajectory 4 * 5 . This is referred to as tangential thrust, because the 
thrusting acceleration will be tangent to the trajectory at all times. 
This methodology provides the maximum increase in the semi-major 
axis of the spacecraft’s orbit over a specific time interval of thrusting. 
Tangential thrust is the nominal thrusting approach used in the 
spiral escape from the departure planet in this program due to its 
near-optimal propellant usage. 

The program allows the user to specify the spacecraft in an 
orbit about the Earth at any date. To generate a trajectory from the 
spacecraft's initial orbit to the final desired orbit about the Moon, the 
plane of the spacecraft's orbit must be aligned with that of the 
Moon s when the spacecraft is escaping from the Earth. This requires 
the spacecraft to transfer between its initial plane of motion to the 
Moon’s plane of motion. Current low-thrust trajectory research has 
not developed adaptive guidance and control algorithms for the 
transfer of a spacecraft between two planes using low-thrust 
propulsion. This lacking instigated the development of a new 
thrusting algorithm for low-thrust transfer between two arbitrary 
planes about the Earth. 

An arbitrary plane of motion in an Earth centered equatorial 
inertial coordinate system can be defined by two of the classical 
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orbital elements, £2 and i, of an orbit in that plane. The longitude of 
ascending node (£2) describes the angle, for a posigrade orbit, from 
the x-z plane to the line where the orbital plane and x-y plane 
intersect. The inclination (/) describes the angle between the orbital 
plane and x-y plane measured at their line of intersection. For any 
orbit about the central body of the Earth, £2 and i will define the 

plane of motion of the object. For two arbitrary planes (see Figure 
8), the angles £2i and i i define the first plane of motion and the 
angles £22 and j '2 define the second plane of motion. The common 
angle between the two planes is called the wedge angle, i'. To enable 
transfers between two arbitrary orbital planes. Plane 1 and Plane 2, 
the angles £2i and i 1 of the spacecraft's orbit. Plane 1, must be driven 
to the angles £22 and 1*2 of the target orbit. Plane 2. This aligns Plane 
1 with Plane 2. In order to change the plane of motion of a 
spacecraft it is necessary to determine how low-thrust accelerations 
would effect the the angles £2 and i that describe the plane. 

The influence of perturbing forces upon the classical orbital 
elements is known from Lagrange's Planetary equations 6 . These 
equations address the rates of change of the classical orbital 
elements due to perturbing forces on the body in orbit. The 
equations that determine effect of perturbing forces on the rate of 
change of i and £2 are: 

iii = F cos ( v + jg) 

dt “ V {8} 

and 

d£2 _ p sin(v + co) 

dt V sin(i) , {9} 

where F n is the normal component of the perturbing force, v is the 
true anomaly of the orbit, co is the longitude of periapsis of the orbit 
and V is the magnitude of the spacecraft's velocity. From equations 
{8} and {9} it can be seen that changes in i and £2 are only caused by 
a perturbing force acting normal to the spacecrafts orbital plane. 

A novel way to determine the necessary thrusting to change 
the spacecraft's orbit from one plane to another was developed using 
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Figure 8 - Two Arbitrary Orbital Planes 
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Lagrange's Planetary equations as a starting point. First, it is 
necessary to develop a desired d//dt and d£2/dt. The desired average 
rates, di/dt and d£2/dt, can be determined by calculating the 
difference (A i) between i\ and 12 and the difference (A £2) between £2 1 
and £2 2 and developing a time period (At) during which the change 
should occur. This At was determined empirically as a function of 
the initial acceleration magnitude of the spacecraft. The At is the 
minimum time required to achieve escape velocity for a spacecraft, 
from starting orbit, under constant tangential thrust. This spiral 
escape was calculated in a nonperturbed two body formulation of the 
equations of motion (see Appendix B). Then di/dt and d£2/dt can be 
approximated by Ai/At and AW/At. The normal force required to 
change i and £2 can be derived from {8} and {9} by rearranging the 
equations to find, 

p . — Ai V 

Atcos(v + co) {10} 


F nn - A£2 v sin(0 
At sin(v + co) , 


where F m - is the normal force to change i 1 to 12 , and F n n is the normal 
force to change £2i to £ 22 . 

Since both £2 and i are affected by accelerations normal to the 
plane of motion, it is not possible for a low-thrust spacecraft to 
change £2 and i independent of each other. Equations {10} and {11} 
are coupled, so the normal forces derived may be in opposing 
directions and effectively cancel during portions of the orbit for a 
given Ai/At and A£2/At. Some previous studies have attempted to 
change only one element, either i or £2, while ignoring the effect the 
of the thrust on the other orbital elements 7 . These coupled equations 
imply that the implementation of the control algorithm that address 
only one of the orbital elements makes an incorrect simplification. 
Additionally, determining the appropriate control algorithm from 
{10} and {11} would be extremely difficult. An essential 
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simplification, however, has been identified to reduce the two control 
equations from {10} and {11} to one single control. 

It should be noted that the assumption that the angle between 
two arbitrary orbital planes can be closely represented by the 
difference between the inclinations, i\ and 1 * 2 , is incorrect. The angle 
between the two planes, commonly called the wedge angle (i*)» is a 
function of both and i. The wedge angle is calculated as follows: 
i’ = cos' 1 [cos(Qi - Q 2 )sin(i'i)sin(j'2) + cos(<i)cos(i 2 )]. {12} 

The control developed for use in this program is based on driving the 
wedge angle between the two planes to zero. The concepts behind 
this control are outlined in the following paragraphs. 

The target orbital plane, in the case of an Earth to Moon 
trajectory, is the Moon’s plane of motion about the Earth. This plane 
is defined as the new x-y plane for the geocentric coordinate system. 
The x-axis points to the Moon at the time of the spacecraft's escape 
from the Earth. The z-axis of this coordinate system is coincident 
with the angular momentum vector of the Earth-Moon system. The 
classical orbital elements of the spacecraft's orbit are then expressed 
with respect to this new frame of reference. Since the desired target 
plane, the Moon's orbital plane, is the x-y plane in the new 
coordinate system, the spacecraft's inclination in this new system is 
the actual angle between the planes or wedge angle. If a thrusting 
control is used to drive the new inclination, or wedge angle, i', to 
zero, then the Q of the spacecraft will also change. However, when i' 
is zero, Q becomes undefined. This is because for any orbit in the x-y 
plane the longitude of ascending node, Q, is undefined. Thus, the 
only control necessary to bring the two planes together is the control 
forcing the inclination, or wedge angle, (i'), to zero (Figure 9). This 
control is identical to equation {10}, 

F n i' = — * 

At cos(v + co) , {13} 

where F n j* is the normal force required to drive the wedge angle to 
zero, v is the true anomaly of the orbit in the new reference plane, 
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Spacecraft’s Orbital 



Figure 9 - New Coordinate System for the Simplified Control 
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and co is the longitude of periapsis of the orbit in the new reference 
plane. 

This control poses practical problems during the integration of 
the spacecraft's orbit. If the inclination of the spacecraft’s orbit is 
zero, or close to zero, integration errors can occur. For this reason the 
integration is performed using equinoctial elements. Atmospheric 
drag is included in the perturbations to the spacecraft's orbit as is 
the oblateness effects of the Earth. These perturbations are added to 
that of the Moon for the entire spiral escape and orbital plane 
alignment. 

When the spacecraft has achieved the required plane change 
the problem is essentially reduced to a planar problem. This is 
because the spacecraft's plane of motion is now the Earth-Moon 
plane and in the absence of out-of-plane perturbations the spacecraft 
will stay in the Earth-Moon plane. The planar problem of cislunar 
transfer has been previously address to some degree 8 . The previous 
work has been extended and modified to account for the three- 
dimensionality of the current approach and the modification of the 
coordinate system to non-dimensionalized coordinates. The control 
of the trajectory generation is passed from the equinoctial 
integration subroutine, EQUIN, to the R3BGEN (Restricted 3-Body 
GENeration) subroutine which integrates the trajectory in the 
restricted three-body problem. 

3.2.2 Midcourse Targeting 

The midcourse targeting portion of the trajectory generation is 
concerned with controlling the spacecraft in order to achieve Earth- 
Moon transfer. The value of the Jacobian constant of a spacecraft is 
used as an indicator of sufficient energy for the cislunar transfer. 

Figure 10 shows a series of zero velocity curves at various 
values of the Jacobian constant. Szebehely notes that when a 
spacecraft in the restricted three-body system has a Jacobian 
constant of approximately 3.3, an equipotential curve like that shown 
at the L -2 point in Figure 10 occurs 9 . When the spacecraft has an 
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energy level that indicates an associated Jacobian constant of less 
than this value, the range of motion of the spacecraft is no longer 
restricted only to orbits about the Earth or Moon, but can include 
transfer orbits between the neighborhoods of the Earth and Moon. 

The midcourse control is based on the desire for the spacecraft 
to have slightly under 3.3 for the value of its Jacobian constant when 
the spacecraft is in the vicinity of the opening in the zero velocity 
curve about the Earth and Moon. This opening occurs at the L 2 point 
shown in figure 10. The spacecraft needs to have the appropriate 
velocity vector in order to pass through the opening between the 
Earth and Moon. It is possible, even likely, that the spacecraft will 
have enough energy for the opening between the Earth and Moon to 
occur, but be in the incorrect position in its orbit to achieve cislunar 
transfer. A control is required that will enable the spacecraft to be 
travelling in the correct direction to achieve Earth-Moon transfer 
when the zero velocity curve opens. 

The Jacobian constant of the spacecraft is calculated at each 
integration step as the spacecraft nears escape from the initial orbit. 
As the value of the spacecraft's Jacobian constant decreases, the 
spacecraft’s energy in the three-body system (potential, kinetic, and 
angular) increases. Figure 11 shows the midcourse control concept. 

It is desired that the spacecraft have enough energy at position 1 so 
that when the spacecraft reaches position 2 the zero velocity curve 
will be open. Different spacecraft acceleration levels correspond to 
different rates of change of the spacecraft's Jacobian constant. This 
means the spacecraft must be very close to escape energy at position 
1 in order to achieve the Jacobian constant of less than 3.3 when the 
spacecraft reaches position 2. The midcourse control is based on 
determining the spacecraft's energy level or Jacobian constant at 
position 1 in order to achieve position 2 using tangential thrust. 

Before the midcourse targeting portion of the trajectory is run, 
the program presents the user with the option of using a parametric 
default value for the Jacobian constant at position 1 or entering a 
different number. The default values of the Jacobian constant are 



Figure 10 - Midcourse Targeting Zero Velocity Curves 




Figure 1 1 - Midcourse Guidance using the Velocity Curves 


parametrically derived as a function of the spacecraft's acceleration. 
This phase of the trajectory generation is user iterative. If the 
spacecraft gains too much energy, it will escape the Earth-Moon 
system. The user is prompted for a response, after 87 days of 
integration time has passed, or if the integration step size falls below 
the nominal value. The program gives the user the options of 
continuing the trajectory generation or regenerating it with a new 
control after Earth escape when either of two conditions is met. On 
trajectories from the Moon to the Earth, the same methodology is 
used. 

The necessary energy for the spacecraft to obtain would ideally 
be the value corresponding to the first zero velocity curve that 
permits travel between the Earth and Moon. In practice, however, 
this is not always the case. The progression of zero velocity curves is 
shown in figure 12. The zero velocity curves start about the Earth 
and Moon as near circular boundaries, Ci. As the energy increases 
the two curves meet at a point in space known as the first Lagrange 
point, curve C 2 . At this position, L 2 , if a spacecraft is placed there 
with zero velocity, with respect to the restricted three-body rotating 
coordinate system, and there are no perturbations, it will remain 
without need for station-keeping. As the energy level increases 
further the two curves join to become a single curve, C 3 , surrounding 
both the Earth and Moon. The shape of this curve is similar to the 
outline of a figure eight. If the energy level increases slightly 
further, the curve, C 4 , surrounding the Earth-Moon system forms an 
opening behind the Moon at Li. This opening means that a 
spacecraft with the appropriate value of the Jacobian constant could 
leave the Earth-Moon system. 

The shape of the zero velocity curves, and the corresponding 
values of the Jacobian constant require consideration when 
developing guidance algorithms in Earth-Moon space. The difference 
between the value of the Jacobian constant where the curves first 
join together and the value of the Jacobian constant when the curve 
permits escape from the Earth-Moon system is not very large. This 


34 




Figure 12 - Zero Velocity Curves in the Earth-Moon System 
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implies that for Earth to Moon trajectories it would be quite possible 
for the spacecraft to obtain a value of the Jacobian constant that 
permits escape from the Earth-Moon system and have the 
appropriate velocity vector to pass through the opening in the zero 
velocity curve behind the Moon. The sensitivity of the trajectory 
generation to the midcourse control value (Jacobian constant) during 
Earth to Moon trajectories is a direct result of this phenomena. For 
Earth to Moon trajectories the difference between the Jacobian 
constant for cislunar transfers and the Jacobian constant 
corresponding to an open zero velocity curve behind the Earth is 
relatively large. Therefore, the sensitivity of the Moon to Earth 
trajectory to the midcourse control value is relatively low. 

3.3.3 Capture and Circularization 

The subroutine CAPTURE controls the capture and 
circularization algorithms during the trajectory generation. The 
capture algorithm is given preference until the two-body orbital 
energy describing a Keplerian orbit about the target planet is 
negative. Then the circularization algorithm is used to lower the 
eccentricity of the orbit. After the eccentricity is less that 0.1, the 
capture algorithm takes over to continue lowering the orbit. If the 
eccentricity of the orbit exceeds the 0.1 limit of circularization 
control again takes control to lower the eccentricity. The two 
algorithms trade control as necessary until the desired final orbit is 
reached. 

As the vehicle approaches the Moon the capture guidance 
phase of the trajectory is initiated. In the absence of impulsive 
thrust, the approach to and capture by the target body are critical 
and must not necessitate maneuvering beyond the limited 
capabilities of the propulsion system. The problem of low-thrust 
spacecraft guidance and trajectory determination between the Earth 
and the Moon was addressed in a study by Richard H. Battin and 
James S. Miller in the late 1950's and early 1960’s. The concept for 
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the spacecraft guidance during capture used in this study is derived 
from Battin and Miller's work 10 . 

The guidance scheme is relatively simple and straightforward. 
The operation of the capture phase guidance is illustrated in figure 
13. The velocity of the spacecraft relative to the target body (i.e., the 
Earth or the Moon) is compared with a parametric velocity profile for 
a reference spiral capture. The parametric velocity profile is a 
function of the radial distance from the target body and the 
magnitude of the thrust acceleration. The difference between the 
spacecraft's velocity and the parametric velocity profile is used in 
combination with the nominal acceleration to determine the direction 
and magnitude of the spacecraft’s thrust vector during capture. 

In order to calculate the parametric velocity profile as a 
function of the radial distance from the capturing body the desired 
reference trajectory must be generated. The desired reference 
trajectory is a spiral trajectory that starts at escape conditions 
relative to the target and achieves circular orbit at the desired final 
altitude about the target body. The subroutine SPIRAL performs the 
generation of the reference trajectory and the calculation of the 
parametric velocity profiles. For a detailed discussion of the problem 
formulation used to generate the spiral trajectory see Appendix B. 

To determine this reference spiral path and the velocity 
vectors that accompany it, the spacecraft starts in a circular orbit at 
the desired final altitude about the target body. The mass of the 
spacecraft at the final altitude is determined by estimating the final 
mass of the spacecraft at the completion of the mission. The 
spacecraft's trajectory is integrated using negative mass flow from 
the propulsion system to spiral out from the target planet using 
tangentially directed thrust. The mass of the spacecraft increases 
thereby decreasing the acceleration as the integration progresses. 
During the generation of the reference trajectory only the 
gravitational field of the target planet is considered. The spiral 
trajectory is otherwise without perturbations and consequently 
remains planar. The calculation of the trajectory continues until the 
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Figure 13 - Capture Phase Thrusting Guidance 
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vehicle is on a parabolic trajectory. The associated radial and 
tangential components of spacecraft’s velocity are recorded at steps 
during the integration as functions of the radial distance from the 
central body. The velocity components are functions of the radial 
distance from the target body. They are obtained by fitting the 
recorded velocity components to polynomial and power curves. 

Figure 14 and 15 are example graphs of the velocity profiles for 
tangential and radial velocity at the radial distance. 

A simplified derivation of the thrust guidance control algorithm 
is presented as follows. This control algorithm determines the 
necessary acceleration vector required for the spacecraft to match a 
desired reference trajectory. The actual velocity of the spacecraft, 

V v , at a given radial distance, r, is compared with the parameterized 
reference capture velocity, V c , at r. The difference between these 
velocity vectors is then determined as V d , where 

V d =V v -V c , {14} 

For a spacecraft flying on the reference trajectory, the 
incremental change in the velocity vector over time, AV C , can be 
approximated as the sum of the spacecraft’s acceleration vector due 
to its thrust and the gravitational attraction of the planet acting over 
some small time increment. At. This implies 

AV c = (a c + g) At {15} 

where a c is nominal acceleration vector of the spacecraft, and g is 
the gravitational acceleration vector of the capture planet. 

The actual change in the velocity vector of the spacecraft, AV v , 
can also be approximated similarly as 

AV v = (a t + g) At {16} 

where a t is the controllable thrust acceleration vector of the 
spacecraft. 

It can be easily seen from equation {14} that the incremental 
change in the difference between the parametric and actual velocity 
vectors, AVd, is simply 

AV d = AV v - AV C . {17} 
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Figure 15 - Radial Velocity as a function of Radial Distance 
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Combining equations {15}, {16}, and {17} and rearranging to solve for 
at, the controllable thrust acceleration vector, equation {18} is 
obtained, 

a t = a c - 

At • {18} 

The thrust acceleration is then chosen so that the rate of 
change of the velocity vector V d is proportional to V d itself. This 
results in 

AV d = V d 

At T c {19} 

where T c is an empirically determined time constant. With this 
formula the appropriate thrust acceleration can be determined in 
both magnitude and direction simply with the knowledge of the 
vehicle’s position, velocity, and nominal thrust acceleration a c . 

In the application of the guidance it is reasonable to assume 
that the direction of the thrust acceleration can be varied at will, but 
the magnitude of the thrust is limited by the capabilities of the 
propulsion system. The spacecraft thrust is never required to 
deliver greater than the nominal thrust. The possibility of a 
reduction in the thrusting acceleration is not precluded as a desirable 
effect of the thrusting algorithm. Figure 16 is a graphical 
representation of the acceleration vectors at and a c . The radii of the 
circles are determined by the nominal acceleration of the spacecraft. 
Then equation {18} becomes, 

a t < a c +-^- 
T c • 

When the appropriate magnitude of the thrust algorithm is greater 
than the nominal capabilities of the engine, a less than nominal 
thrust is used in the appropriate direction. 

After the spacecraft has been captured by the target body, the 
capture algorithm continues to lower the spacecraft's energy level to 
bring it in closer to the desired orbit. One distinct difficulty of the 
capture algorithm is lowering the eccentricity of the capture orbit. 
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For this reason a separate circularization algorithm was developed to 
control and lower the spacecraft’s eccentricity during capture. 

The concept behind the circularization algorithm is very 
straightforward. The spacecraft needs to be in a near circular orbit 
for the capture algorithm to be most effective. When a orbit is 
nearly circular the velocity at each portion of the orbit is 
approximately the same magnitude. For eccentric orbits the velocity 
of the spacecraft is highest at periapsis and lowest at apoapsis. The 
goal of the algorithm then is to drive the spacecraft’s velocity to a 
uniform value along the entire orbit. At periapsis the spacecraft 
needs to increase its velocity, at apoapsis it needs to decrease its 
velocity. In order to achieve this, the spacecraft’s thrusting vector is 
pointed opposite its velocity vector at apoapsis and along the velocity 
vector at periapsis. This thrusting behavior translates into the 
following controls for the radial a tr , and transverse a ts , acceleration 
components. Then 

a tr = a c sin(v) 

and 

a, s = a c cos(v) 

where v is the true anomaly of the orbit. Figure 17 shows what the 
acceleration vector would look like at various points along the orbit. 

A potential restriction on nuclear-powered OTVs is the 
proposed nuclear safe orbit (NSO) 11 . This would be a designated 
altitude above the Earth, below which the nuclear powered 
spacecraft would be prohibited. The spacecraft would be unable to 
descend below the restricted altitude at any point of its trajectory. 
This limits the types of trajectories available for cislunar transfer. 

1 Advanced Space Analysis Office-Sverdrup/NASA-LERC, “Evaluation of 
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Appendix A 


Variation of Parameters 

This is a useful method when a solution is sought to a 
homogeneous differential equation under a perturbation. The 
homogeneous equation under consideration is the equation of motion 
of a spacecraft about the Earth. The equation of motion is commonly 
written as, 

- _ px 

' x 3 « {Al} 

where p is the gravitational parameter of the Earth, x is the position 
vector of the spacecraft, and x is the resultant acceleration of the 
spacecraft. This equation has a homogeneous solution that derives 
the six classical orbital elements a a = (a, e, i, M, to, Q). An alternate, 
but equally valid solution derives the equinoctial elements p a = (a, h, 
k, A,, p, q). This means the position and velocity of the spacecraft is a 

function of time only, so x = x(t), and x = x(t). 

If there is a perturbing force, F, to this orbit, the equation of 
motion can be written as 


** px , r: 

x 3 . {A2} 

this is an associated inhomogeneous differential equation and the 
solution to equation {A2} can be found through variation of 
parameters. 

The basic idea behind variation of parameters is to replace the 
constants of the homogeneous solution (the orbital elements) with 
time-varying functions. This implies the position is a function of 
both the orbital elements and time, x = x(p a ,t). Then to determine x, 
the chain rule for differentiating is used. First the partial of x with 
respect to p a is multiplied by the partial of p a with respect to time, 

and then the partial or x with respect to time is taken, so 

^ _ 3x 9p a dx 


dp a 3t 9t. 
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However, this equation is constrained by setting 

dx a P« _ Q 

dPa , {A3} 

this forces the velocity to be the same for the perturbed and 

unperturbed case. In effect, this constraint says that the position of 

the orbit is not changed at the time the perturbation acts. To 

determine the acceleration, x, the chain rule is again used, 

^ 3x 3p a dx 
x = — + — 

9p a dt 9t. 


However, equation {Al} already gives a solution to x. This implies, 

ax 8p a _ p 

a Pa at . {A4} 

The six unknowns of equations {A3} and {A4}, p 0 , can be found 


by inverting the 6 x 6 matrix of partials. If equations {A3} and {A4} 
are put in a matrix, the inverted 6x6 matrix of equations would 
yield. 




. [o 

dx 

dx 

LF- 


Pa = 


If the matrix equation is multiplied through, it is found, 

p«=%-F, 

9x 


{A5} 


where the right hand side of { A5 } is the dot product of the partials 
with respect to velocity together with the rectangular components of 
the perturbing acceleration F. To solve for p« it is first necessary to 
know the explicit form of the partial derivatives of the orbital 


elements with respect to the velocity components of the reference 
orbit. These partial derivatives are obtained by using the Poisson 
Brackets 1 : 


dpq 

dx 


• £ (Pc Pf,) §- . 

p = 1 dpp 


{ A6} 
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To use the Poisson Brackets to find the partials of the 
equinoctial elements with respect to the velocity, the partial 
derivatives of position and velocity vectors with respect to the 
equinoctial elements must first be found. The initial step to this is to 
determine the position and velocity of the spacecraft from the 
equinoctial elements of the orbit. The quantities (Xj, Yi, 0) are the 
coordinates of the spacecraft relative to the equinoctial frame. These 
coordinates must be found and then transformed into the traditional 
(x, y, z) of the inertial coordinate system. 

It is necessary to recall the coordinate system for equinoctial 
elements. The equinoctial coordinate frame is defined by unit 

A XV. ^ 

vectors f, g, and w illustrated in Figure A1 and defined by 

1 - p 2 + q 2 

f = \ 2 2 2pq 

-2p J, 

2pq 

1 + p 2 - q 2 

2 q 

2p 
-2q 

i - p 2 - q 2 . . 

In order to have a variation of parameters program which is 
valid for all eccentricities and inclinations it is necessary to use a 
formulation which is free of singularities. The key to this is to have 
the angles describing the orbit defined for all eccentricities and 
inclinations. The equinoctial formulation uses angles called the mean 
longitude (X), and the eccentric longitude (F), rather than the more 
classical angles of mean anomaly (M), eccentric anomaly (E), and true 
anomaly (v): 

X = M + co + Q, { A7 } 
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Figure A1 - Equinoctial Coordinate Frame 
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F = E + cd + £2, {A8} 

The position in the orbit can be indicated by F. However, in order to 
compute the position vector of the spacecraft from the equinoctial 
elements, it is necessary to solve Kepler's equation. It is 
advantageous to write Kepler’s equation in terms of the eccentric 
longitude rather than the eccentric anomaly because the eccentric 
longitude is defined for all eccentricities and inclinations. Some 
elementary manipulations show the Kepler's equation and the 
expression for the radial distance, r, from the Earth can be written in 
terms of the eccentric longitude, 

X = F + h cos F = k sin F, { A9} 

r = a [1 - h sin F - k cos F]. {A10} 

Kepler's equation {A9} can be solved for F with the standard 
Newton-Raphson procedure (or any other iteration method) once the 
value of X has been determined from {A7}. Once Kepler’s equation 
has been solved, the three coordinates (x, y, z) of the inertial system 
are obtained with the following matrix equation: 


x 

y 

z 


1 

1 + p 2 + q 2 


1 


p 2 + q 2 



' X! ' 


Yi 


. 0 . 


2pq 2 p 

2pq l+p 2 -q 2 -2q 

-2p 2q l-p 2 -q 2 J 

The quantities (Xi, Yi, 0) are the coordinates of the spacecraft 
relative to the equinoctial frame. They can be expressed in terms of 
F by: 

h (X - F) 


Xi = a 


cos F - k + 


i + vr 


+ e J 


Yi = a 


sin F - h + 


k (X - F) 

1 + VT 


+ e 2 J. 


Also of note are the following time derivations: 

F = HA 
r ’ 


{All} 





x, = ^oa 

r 


2 


. tj h (k cos F + h sin F) 

(1 + Vl - e 2 ) 



cos F - k ( k cos F + h sin F) 

(1 + Vl - e 2 ) 


The above expressions are valid for all eccentricities and 
inclinations. In addition, all the expressions are functions of the 
equinoctial rather than the classical elements. Due to the lack of 
singularities for circular or zero inclination orbits the equinoctial 
formulation is more appropriate for the integration of low-thrust 
trajectories. 


The partial derivatives 

In the variation of parameters derivation it is necessary to 
have several partial derivatives of the two-body equations. First, all 
of the partials of Xi and Yi with respect to h and k are needed. The 
following expressions were taken from work by Broucke 2 : 


axi 

3h 


= a 


- (X - F) (p + 


-22_) - 3- cos F (hp - sin F) 
1 - p r 


axi 

3k 


= -a 


(X-F) 


—^2— + 1 + 3- sin F (sin F - hB) 
1 - p r 


?Yi 

3h 


= a 


(X-F) 


^2- -1+3. cos F (kp - cos F) 
1 - p r 


and 


3Yi 

3k 

The quantity P 


ST 

• 

1 

k 2 B 3 
P+ P 

+ 3- sin F (cos F - kp) 

L 

1 - PJ 

r 


which has been introduced here is defined as follows: 

P= 7 1 

1 + V 1 - h 2 - k 2 • 
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With the use of the above expressions the following partial 
derivatives of the position vector x are easily derived: 

— = -J- [x - 2_x (t - to)] 

8a a 2 

8x 8Xi z 8Xi - 

— = — - f + — L g 

8h 8h 8h , 

8x 8Xi - 8Yj - 

— = — l - f + — — g 

8k 8k 8k , 

8x _ Xi f+ Yi g 

8X 0 n 

— = \ [q (Yif-XiS-Xj w] 

8p 1 + p 2 + q 2 , and 

— = \ [p (Xig-Y,f)-Y, w] 

8q 1 + p 2 + q 2 

Finally the partials of the equinoctial elements with respect to the 
velocity are obtained by using the Poisson Brackets shown in 
equation {A5}. The results are then: 

da _ 2x 
8x n 2 a , 

8q _ (1 + p 2 + q 2 ) Xi — 

dx 2na 2 Vl - h 2 - k 2 

dp _ (1 + P 2 + q 2 ) Yj — 

dx 2na 2 V 1 - h 2 - k 2 


5 3 



8h 

dx 


k(q Y t - pXQ 
na 2 Vl - h - k 2 


w + Vl - h 2 -k 2 

n q 2 



h pXi)? + (§Xi 

n / Uk 



g 


ak 

dx 


-h(qY! - pXQ - 
na 2 Vl - h 2 - k 2 




- k pXi)f%(3Yi 

n I \ah 


+ kB^i- 
n 


g 


3^0 _ -2 


9x 


na- 


x 2-x (t - to) 




na- 



ll 1 

Ll ah 

+ k WM 
ak j 


i 

[c 


fH' 


gYi 

ah 


+ k 


aYi 

a k 


These partials are calculated at each integration time step and 
used with equation {A5} to find the time rate of change of the 
equinoctial elements due to the perturbing force F. This permits the 
integration of the equinoctial orbital elements under perturbation. 


g 


1 Broucke, R., and Cefola, P. “The Equinoctial Orbit Elements,” JPL Technical 
Memorandum 391-238, September 28, 1971. 

2 Broucke, R„ and Cefola, P. “Numerical Integration of Satellite Orbits with 
Equinoctial Elements,” JPL Technical Memorandum 391-248, November 5, 1971. 
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Appendix B 



The equations of motion of a powered spacecraft in a two- 
dimensional orbit can be derived the spacecraft's kinetic energy, 
potential energy, and from the force equation, F = M*a (total force = 
mass * acceleration). 

If the magnitude of the acceleration of the spacecraft due to 
the propulsion system is a, then 

T = M*a, {Bl} 

where M is the mass of the spacecraft. The spacecraft’s acceleration , 
a, is then determined from equations {Bl} to be 

a=X 

M 

The acceleration of the spacecraft due to the thrust of the 
propulsion system changes as the total spacecraft mass, M, changes. 
As mass is expelled from the spacecraft the total spacecraft mass 
decreases. So, 

M = M 0 - m*t , 

where t is the time period during which the spacecraft is losing mass 
at m, the mass flow rate of the spacecraft's propulsion system. 

To determine the equations of motion of the spacecraft in polar 
coordinates the radial and transverse components of the acceleration 
due to the thrust need to be calculated from the thrusting force. 
Figure Bl shows an arbitrary thrust vector with respect to the radial 
and transverse vector. Assuming that the thrust vector makes an 


angle $ with the local horizon, (measured from the direction of the 
motion), the radial, R, and transverse, S, components of the 
acceleration are: 


a R = -I-sinA, ac = r-^cosd> 
R M S M 


In an unperturbed orbit the only external force present is the 
Earth's gravity force which in polar coordinates has the radial 
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component. The equations of motion in polar coordinates can be 
derived quickly from the Lagrangian which is the sum of the kinetic 
and potential energy, 

L = -L (r 2 + r 2 0 2 ) + . 

2 r 

Using the Euler-Lagrange equations, the equations of motion are 
obtained as, 


r* -rt = -GM 

r-2 


and 


d(r 2 6) _ Q 

dt ' {B3} 

If the acceleration components due to the spacecraft's thrust from 
equation {B2} are added to the equations of motion {B3 } , then 


r - r0 = - GM + X. sin<t> , 


M 


{B4} 


= r i cost, 
dt M 

The above system { B4 } can be replaced by a new system of 
four first-order equations. This is done by introducing two new 
variables, the radial velocity component u = r, and the transverse 
velocity component, x> = r0. Through the change of variables, {B4} can 
be rewritten as, 

u - u0 = - OM j- T 

T 4 

{B5} 

= rX- cos<|) 
dt M 

The two second order differential equations, {B5}, can be broken into 
four first-order coupled equations: 

r = u , 


+ — sin<J> , 
.2 M 


e-?- 
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u = ^ + X s in(j) 

r r 2 M 

i) = -W-uu + iX cos6) . 
r M 

For the numerical integration of these four equations, the state 
vector is defined to be (r, 0, u, u), (in this order). These equations 
were used in the subroutine SPIRAL to perform generation of the 
velocity component profiles for the reference trajectory for spiral 
capture. This formulation was especially useful because the 
tangential and radial components of the velocity were required to be 
save to an array as a function of the radial distance from the planet. 
Using this formulation eliminates the need for any cumbersome 
calculation of the velocity components and radial distance as the 
integration progressed. The necessary information was inherent in 
the formulation of the problem. 
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